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l-known methods for solving the shape-from-shading problem require 
t reflectance map. Here we show how the shape-from-shading problem 
ten the reflectance map is not available, but is known to have a given 
unknown parameters. This happens, for example, when the surface is 
libertian, but the direction to the light source is not known. We give 
rithm that alternately estimates the surface shape and the light source 
f the unit normal in parameterizing the reflectance map, rather than 
stereographic coordinates, simplifies the analysis. Our approach also 
flve scheme for computing shape from shading that adjusts the current 
ocal normals toward or away from the direction of the light source. The 
tment is proportional to the current difference between the predicted 
l brightness. We also develop generiilizations to less constrained forms 
ps. 
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E, and a reflectance map, R, the shape-from-shading problem may be 
of recovering a smooth surface, z(x,y ), satisfying the image irradiance 


E(x,y) 


dz dz 
dx ’ dy 


1, of E. Any given boundary conditions on z{x ) y) should also be satisfied. 
;es the form of a first-order partial differential equation (Horn, 1970 & 


is formulation are a number of assumptions, the principal one being that 
‘ a surface patch does not depend on its position in space. Another is 
spicts a smooth surface of homogeneous reflectance. Several algorithms 
1 to tackle the problem, notably those of Horn (1975), Strat (1979), and 
1981). 

my difficulties these schemes face in practice is that the reflectance map 
lown. The reflectance map specifies how the brightness of a surface patch 
dentation under given circumstances. It therefore encodes information 
ing properties of the surface, and information about the distribution 
the light sources. In fact, the reflectance map can be computed from 
reflectance-distribution function and the light source arrangement, as 
i Sjoberg (1979). 

tering a new scene, the information required to determine the reflectance 
>t available. Yet without this information, the shape-from-shading prob- 
mulatcd, much less solved. The dilemma may be resolved if a calibration 
shape appears in the scene, since the reflectance map can be computed 
Iere we wish to consider the situation where we are not that fortunate, 
ng to evaluate how some basic assumptions can resolve this impasse, 
has looked at the problem of recovering shape from shading under the 
the image depicts a Lambertian surface illuminated by a point source 
3 unknown. Under the additional assumption that the surface is locally 
: normals are shown to be recoverable by a local operation. This method 
on the iterative propagation of information across the image, 
e serious drawbacks to the load approach, however. One problem is that 
imption is very restrictive. In fact, spheres are the only surfaces whose 
bilical. So this method naturally computes incorrect normals for other 
the errors for approximately spherical surfaces, such as ellipsoids of 
nay be acceptable. Further, the constraining effect of known occluding 
Is cannot be incorporated into the local method. This is unfortunate 
finals provide powerful boundary conditions on the shape-from-shading 
n by Bruss (1983). Finally, because a local method does not take into 
s when calculating the normal at a point, nearby normals may differ a 
ularly in the presence of noise. 

it an alternative approach that does not suffer from these disadvantages. 
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! scheme for shape and source direction 

jeover the shape of a smooth surface depicted in an image, E , that is 
gion fi in the zy-plane. Let the shape of the surface be characterized 
n, that associates a unit normal with each point in 0. The problem 
[id n(x,y) over 0. Assume for now that the object has a Lambertian 
it is illuminated by a single point source. If the vector s points to the 
[/) is the unit normal of a surface patch, then the apparent brightness of 
[i by the reflectance map 

R a (n(x,y)) = n(x,y) ■ s. 

e way, force s to be a unit vector; this allows for the possibility that the 
pirce may be unknown. This way we can also deal with unknown sensor 
liknown surface albedo, provided it is uniform. 

then becomes one of finding a smooth shape, n, and source direction, s, 
fige irradiance equation 

E(x,y) = n(x,y) • s V(x,y) E fh 

rightness cannot be determined with perfect accuracy, and so it appears 
nsforrn this into a minimization problem (Ikeuchi & Horn, 1981; Horn 
There is another reason to consider this as a minimization problem: 
to solve the image irradiance equation as it stands, we obtain a set of 
ions equivalent to the characteristic strip equations. Here, however, we 
scheme lending itself to a parallel implementation on a grid, as originally 
n (1970). Further, a shape-from-shading problem that has noisy image 
•t have a theoretical solution. A minimization approach will, however, 
ry of a shape that fits the given data best, in a sense determined by the 

Doth shape, n, and a source direction, s, that minimize 

Jj (E{x,y) - n(x,y) • s) 2 dxdy. 

is, and there are no errors in brightness measurements, then the image 
on will have been satisfied (although there is no guarantee that the 
jc integrable; see Horn & Brooks, 1985). 

a regularizing component (Poggio & Torre, 1984) by incorporating the 



to select a particularly smooth solution from a possibly infinite set of 
that a subscript here denotes partial differentiation, and that squaring 
cut to taking the dot-product with itself. Finally, we wish to insist that 
length. This is accomplished with the constraint 


n 2 (x, y) 1 v(z,t/)efi. 
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iree terms gives the composite functional 

) = ff (£-n-s) 2 + A(n“-f n 2 y ) +/i(x,y)(n 2 - 1) dxdy 

inimized with respect to n and s. Here, A is a scalar that weights the 
ice of the regularization term, while n(x,ij) is a Lagrangian multiplier 
impose the constraint that n (x,y) be a unit vector (see Horn k Brooks, 

is a problem in the calculus of variations. First, assume that s is known 
be minimized by a suitable choice of n. Extrema of functionals are 
the associated Euler equations (see Courant and Hilbert, 1953). The 


ll F(x,y,n,n Xy n y ) dxdy 


ation 


d d 

p _ p _ p — n 

1 n 0 r n s ~ r n — u. 

ox ay 

m, it follows that I has the Euler equation 
{E — n • s)s -f AV 2 n — pn = 0, 


o d 2 d 2 

=r-f-- 

dx 2 dy 2 

pproximation to the Laplacian operator is given by 

NMr, ~ ~2^l} - n 0')> 

istance between adjacent picture cells in the image, and rqy is the local 
s 

n ij ~ l( n i',y+.i + n *j-i + n i+i,j + 
mslate the Euler equation into the discrete form 


'tA 

\fiij ~h y — Pij — 0* 

in order to isolate n,y on one side yields the iterative scheme 

= I : J 4\) (”£» + TX^' 3 ~ ' S ^)’ 

ape, given the light source direction. Other approximations for the 
id to improved results, at the cost of increased computation. For ex¬ 
lie more accurate 9-point approximation for the Laplacian, in which 
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'o H + n «,y-j + n.+ij + n »-ij) 

+ + n i-i,y+j + n i+ij+i + n !+ij-i) 


then twice a 5 many array accesses are needed (and the constant multiplier e 2 /4A 
becomes 3e7l0. ). The simple 5-point approximation was adequate for our purposes. 

Note that wc have yet to solve for /i, the Lagrangian multiplier. This can be avoided, 
however, by obsi rving that the division of the right hand side by (l f/i(e 2 /4A)) does not 
change the direction of the vector being computed. Since fi is intended to ensure that 
the result is nor ualized, we can simply do this explicitly, as in 


C - 

in 1 1 in* 1 

IJ / I l] 


n f, ■ s ) s 


Now conside: the problem of minimizing I with respect to s, given that n is known. 
This is a proble ti in conventional calculus. Computing the partial derivative of I with 
respect to si we lave 


2 [E — n • s)n dx dy = 0, 


and so 


E n dx dy + / / (n • s)n dx dy — 0. 


Noting that 


(n • s)n — ( n T s)n = n(n T s) — (nn r )s, 


then by substitution we have 


' Ei 1 dxdy — / / nn dxdy s, 

n IJJn 


where (nn 7 ), aJd also the integral of (nn 7 ), arc 3x3 matrices. From this we finally 
obtain the c esircll equation 


nn dxdy 


E n dx dy. 


Here 1 denotes the inverse of a matrix. A discrete version of this formula, in which the 
integrals are replaced by sums, is easily obtained. An iterative scheme in both n and s 


3. Properties and il'rformanee of the scheme 


is then witliin giasp: 
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This takes advai tage of the fact that s need not he a unit vector. If the brightness of the 
source is kr own we cati normalize s so that it is a unit vector. Then the determination 
of s becomes slq fitly more complex, since it involves a constrained minimization. 


3. Propertlies Ind performance of the scheme 
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leme has two components: one concerned with the recovery of shape, 
led with the determination of the source direction. The shape-recovery 
n intuitively satisfying form. In essence, a new normal is computed by 
erage, and adjusting this either toward or away from the source. The 
gn of the adjustment is determined by the brightness error of the current 

:iape, a new source direction is computed by a single pass through the 
pe-recovery, no iteration is necessary. The 3x3 matrix (mi 7 ) is summed 
as is the vector [E n). The source direction can then be computed using 
tion or even Cramer’s method (sec Korn and Korn, 19G8). The sourcc- 
ent has been tested on a number of images and shapes. When the data 
the estimate of source direction is extremely accurate. Furthermore, 
very good in the face of significant noise. For example, a synthetic image 
a sphere illuminated by a point source in the direction ( — 4, 3, 8) r . The 
ized to 255 irradiance levels, and the correct surface normal was given 
250 image points. Gaussian noise was added to the image giving an 
tion in irradiance values of 34. Despite this, the source-finder computed 
•urce direction that was only 2.7° in error. Further trials gave similar 

rection estimates are robust because the whole image is used. Theo¬ 
rem is highly over-determined, as source direction is recoverable from 
of only three different surface orientations. Using the whole image 
that noise effects arc significantly reduced. 

J-sourcc-from-shading problem for a point light source and Lambertian 
■ural two-way ambiguity. It the image irradiance equation is satisfied 
ipe ill and the source direction si, it will also be satisfied by the dual 
rce direction S 2 where 


ri 2 = 2z(ni • z) - ni 
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Shape and source from shading 

wing direction. Here, both source direction and surface normals are 
he viewing direction. This is easily verified by observing that 

2z(ni • z) - ni) • (2z(si • z) - si) 

(z • z)(ni • z)(si • z) - 2(ni • z)(si • z) - 2(s A • z)(n x • z) + ni • si 

ii si. 

al values for the normals, the shape-and-source scheme will head for 
of these solutions. The dual shape and source direction can then be 
:diately using the equations given above. 

snt two examples of the program at work. Each (synthetic) image con- 
a Lambertian surface illuminated by a point source in the direction 
lages each contained more than 1000 points at which normals were to 
lormals were assigned an initial value of (0,0, l) r , as was the source di- 
ig boundary normals were given. The equations for n tJ could be solved 
l the Gauss-Seidel algorithm. Since we are ultimately interested in par- 
ion on a grid, we used the Jacobi method instead (despite the fact that 
method has slightly better convergence properties), 
ge portrayed a hemisphere viewed from directly overhead. After 100 
= 0.005, the average angular difference between estimated and correct 
than 3°. The maximum such deviation was less than 2.5 times the 
’he estimate of the light source direction, at this time, had errors in 
th angles of 1.4° and 1.6° respectively. 

lage depicted a cylinder with rounded, hemispherical ends, viewed from 
ndicular to its axis. After 60 iterations, this time with A = 0.003, the 
rror in surface normal was less than 5°. A further 30 iterations brought 
) 4°. The maximum error remained somewhat larger, however, due to the 
y to smooth the intersection between the cylinder and the hemispheres, 
mutli and zenith angles for the source were 7.3° and 1.1° respectively, 
iterations. These, too, improved slowly with further processing, 
as sometimes slow in converging. After rapid initial improvements, the 
ould decrease appreciably. However, one might expect the scheme to be 
of the current iterative methods, given the disadvantage of not knowing 
direction. Convergence could be accelerated by employing multigrid 
ropagate information across the image more quickly (see Terzopoulos, 
gly, in the examples considered, a reasonable estimate for the source 
ained after only a few iterations. Subsequent processing just improved 

ernes for other reflectance maps 

wo more new iterative schemes: the first extends the shape-and-source 
uations in which a simple model of the sky is also included; the second 
r eloped above to find shape from shading, given a general reflectance 
t recover source direction. 
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R$ky (n) = \{l + n • z) 

ation in which a Lambertian surface is illuminated by a hemispherical 
i radiance (Brooks, 1978; Horn h Sjoberg, 1979). A point source may 
map to give 

R. iS (n) = a(n • s) + ^(1 + n • z) 

the relative intensity of the sun and the sky. We can now generate a 
-and-source recovery, under the assumption that the image was formed 
th the reflectance map R as . 
inimize 


JJ T-a(n-s)-^y(l + n- z)| + A(n^ + njj) + fi(x, y)(n 2 - 1) dxdy , 
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>oth n and s. Fixing s for the time being, we are required to minimize 
>nal with respect to n alone. The Euler equation for this problem is 

- a(n • s) - ^r(l + n ‘ z)) (as + ^^z) + AV 2 n - /m = 0. 

‘ore, the following scheme is obtained: 
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Thus we o 


ectjnce map, R sa , has been substituted back into the equation to improve 

ue n to be fixed and minimize the functional with respect to s. This we 
differentiating with respect to s and equating the result to zero. Thus 

-2 JJ - o:(n • s) - -lp(l |n-z)Jan dxdy = 0. 

J (F- ^(l + n-z))n = JJ a(n • s)n dxdy. 

e that (n • s)n = (nn r )s, the equation becomes 

J (E — (1 + n • z)) n dxdy = a J J nn T dx dy s. 

e equation in s given by 

j — JJ nn T dxdy^ JJ (E — ? (1 -f n • z)) n dxdy. 
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is required from shape. Indeed, if used in this way, the shape-recovery 
be generalized to incorporate any reflectance map, R(n). In minimizing 
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te in discrete form and combine with the iterative scheme for n derived 
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|ips fie most appealing of the current shape-from-shading schemes that deal 
ectance map. It is simply derived, and is expressed elegantly in terms 
ather than a two-parameter system such as the stereographic fg space 
>rn. 


ve use of the unit normal constraint 


ving the shape-and-source finder, we avoided solving for the Lagrangian 
. Instead, we normalized the estimate for n after each iteration. We 
me in which the multiplier is dealt with explicitly, 
ninimize the functional 
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thejEuler equation 

(E - R)R n + AV 2 n — /m = 0. 

[lot roduct of this with n we find that 

H = (E-R)(R a - n) + A(V 2 n-n), 
itut|ng for f-i in the Euler equation, we get 

{E - R) [R n - (R„ ■ n)n] + A [V 2 n - (V 2 n • n)n] = 0. 

(x • n)n = (n r x)n = n(n r x) = (nn r )x, 

|or x| and letting Mbethe3x3 matrix 

M = I - nn r , 

uatiln reduces to 

M [{E - R)R n + AV 2 n] = 0. 

in components orthogonal to n, since 

(I - nn r )n = n — ( nn T )n — n - n(n r n) = n - n(n • n) = 0. 

s provides only two constraints on the solution vector n. The remaining 
= 1. Note that, because M is singular, we cannot simply eliminate M 
l above by multiplying through by its inverse, 
ard finite difference approximation such as 

V 2 n « ^(n - n), 

k loiil average of n given earlier, we can write the Euler equation in the 
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, we omit subscripts.) We can then develop an iterative scheme in 
|due, m, say, is used for the center term in the above approximation, 
ms are computed using the old value of n. This way we obtain 
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^n, where p J_ n. Then m 2 = p 2 + v 2 and M m = M p = p, since 
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are theoretically two solutions for i/, one positive and one negative. The 
ids to a new estimate close to the previous one, while the negative value 
almost opposite to the old one. It is clear that one should use the positive 

lly have the scheme 

p‘ +1 = (I - ) (n*. + ~(E t] - ))/?„(n‘)) 

=pf, +I +<V‘ - ( p f» +1 ) 2 - 

pled with the source-recovery component given earlier, when /?(n) = n-s. 
tch a solution with this scheme, p will be small, since n ^ n and E & R. 
early stages of iteration, it may be necessary to place an artificial limit 
f adjustment made away from the old normal. That is, one may have to 
ide of p so that problems do not arise in computing u using the equation 

nethod for solving the underdetermined equation Mm = 0, under the 
1, can be arrived at most easily using the pseudoinverse of the matrix 
this approach in the exposition here since the solution can also be found 


thods for obtaining shape from shading assume complete knowledge of 
lap. Here, we considered the situation in which a Lambertian surface 
r a point light source from an unknown direction. Thus we dealt with 
rather than a fixed, reflectance map. The local approach to recovering 
ing, which is also intended to deal with unknown source direction, was 
Jveral drawbacks, notably its restrictive assumption that surfaces are 

of unit normal vectors for describing surface orientation was important 
ent of our method. This led to simple derivations and elegant presenta- 
:m was cast as one of minimizing a positive-definite functional containing 
r term, a regularizing term, and a Lagrangian multiplier to enforce the 
e normal be of unit length. The Euler equation for this calculs of varia- 
shown to be a second-order partial differential equation in the unknown 
metion. A convergent iterative scheme solved it in the discrete domain, 
of the light source can be determined in closed form if the surface shape 
g any iteration a source-direction estimate can be obtained using the 
of the surface shape. The iterations for obtaining increasingly accu- 
the surface shape can be interlaced with estimation of the light-source 
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